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The Thermal Structure of Gravitationally-Darkened Classical Be 

Star Disks 

M. A. McGill\ T. A. A. Sigut\ & C. E. 3ones\ 
ABSTRACT 

The effect of gravitational darkening on models of the thermal structure of 
Be star disks is systematically studied for a wide range of Be star spectral types 
and rotation rates. Gravitational darkening causes a reduction of the stellar ef- 
fective temperature towards the equator and a redirection of energy towards the 
poles. It is an important physical effect in these star-disk systems because the 
photoionizing radiation from the central B star is the main energy source for 
the disk. We have added gravitational darkening to the bedisk code to produce 
circumstellar disk models that include both the variation in the effective tem- 
perature with latitude and the non-spherical shape of the star in the calculation 
of the stellar photoionizing radiation field. The effect of gravitational darkening 
on global measures of disk temperature is generally significant for rotation rates 
above 80% of critical rotation. For example, a BOV model rotating at 95% of 
critical has a density-averaged disk temperature ~ 2500 K cooler than its non- 
rotating counterpart. However, detailed differences in the temperature structure 
of disks surrounding rotating and non-rotating stars reveal a complex pattern of 
heating and cooling. Spherical gravitational darkening, an approximation that 
ignores the changing shape of the star, gives good results for disk temperatures 
for rotation rates less than ^ 80% of critical. However for the highest rotation 
rates, the distortion of the stellar surface caused by rotation becomes important. 
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Introduction 



A classical Be-star is a non-super giant B star that possesses a gaseous equatorial disk 
( iPorter &: Rivinud |2003[ ) . This circumstellar material produces an emission line spectrum 
including a prominent Ha emission line. Some Be stars, for example 7 Cassiopeia, produce 
a rich emission line spectrum with many hydrogen lines and some metal lines. Be star disks 
also produce a continuum excess due to bound-free and free-free emis sion. This excess has 
been observed at a va riety of wavelengths from the visible to the radio (ICote fc Waterslll987l : 
Waters fc Cotg 119871 ). Linear polarization of up to 2% has been observed in the continuum 
emission and is caused by electron scattering in the non-spherical circumstellar envelope 
( iMcLean fc Brown! Il978t iPoeckert et al.lll979l ). Be star disks have b een partially resolved 
by lo ng-baseline interferometric measurements, initi ally in the radio (iDoughterty fc Tavlor 
I992I ). and later at infrared and optical wavelengths (IStee fc Bittarll200ll : iTycnerl 120041 ). 



The stellar absorption lines of Be stars indicate rapid rotation of the central star, 
and t his property is thought to play an important role in the release of material into the 



disk (IPorter 



19961). There is strong evidence that Be star disks are in Keplerian rotation 



( iHummel fc Vranckerul2000l : lOudmaijer et al.ll2008l ): however, it is still uncertain if the equa- 
torial rotation speeds of Be stars ar e close to "crit ical" , i.e. close to the Keplerian orbital 
speed at the inner edge of the disk ( ICranmeii l2005l ) . In addition, another mechanism must 
somehow allow for the wide range of behaviours observed in Be stars: some systems have 
exhibited essentially stable Ha emission for as long as they have been obse rved, while others 
are h ighly variable with timescales ranging from less than a day to decades ( IPorter &: Rivinus 
2003h . 



A rapidly rotating star is changed in two ways: the stellar surface is distorted with 
the equatorial radius becoming larger than the polar radius, and the stellar surface tem- 
perature acquires a dependence on stellar latitude with cooler gas at the equator com- 
pared to the hotter pole. These phenomena together are commonly called gravitational 
darkening. The rotational distortion and surface intensity variation with latitude have di- 
rect observational support from interferomet ric observations, incl uding the stars Achernar 
( iDomiciano de Souza et al.ll2003l ) and Altair (IMonnier et al.l 120071 ). 



Because the equatorial region of the photosphere is both the fastest rotating and the 
dimmest, it can be difficult to detect this maximally Doppler-shifted region in integrated 
light. This problem can be compounded by the obscuration of the equatorial regions 
by disk material. As a result, the bounds on Be star rotation rates remain contentious 
(IChauville et al.lboOli rTownsend etaP 120041 : [Cranmei]l2005l : IPremat etauboosl ). Neverthe- 
less, the canonical result is that most Be stars rotate at ~ 80% of their critical speed (as 
defined in Section [2]). Unfortunately, it seems that directly measuring projected rotational 
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velocities above 80 % of the critical speed may not b e possi ble due to the dimming of the 



equatorial regions ( iTownsend et al.l 12004 ) . ICranmerl (|2005[ ) presents a detailed statistical 



analysis of Be star v sin i values in which he attempts to match the observed distributions 
using a parametrized distribution of equatorial velocities. He concludes that early-type Be 
stars, 07-B3.5, must rotate at rates significantly less then their critical speeds (peaking 
at 40-60%) while later type Be stars, B3.5-A0, rotate much closer to their critical speeds 
(peaking at 70-90%). 

Previously, the effects of gravit ational darkeiiing ha ve been includ ed in disk models fo r 
the stars Achernar, or a Eridani, ( ICarciofi et al.l l2008l ) and ( Tauri ( ICarciofi et al.l 120091 ). 
However, the rotating star was assumed to be a spheroid rather than following a Roche 
profile (see Section |2]). In the literature, gravitational darkening has been included only as a 
fit parameter in models of individual stars. A comprehensive study on the how the inclusion 
of gravitational darkening affects models of Be star disks has not yet been performed. 

In this paper, we systematically examine the effect of gravitational darkening on the 
thermal structure of a circumstellar disk. We restrict our investigation to changes in the 
energy reaching the disk and the resulting changes in the temperature distribution. Sec- 
tion [2] outlines the basic theory behind gravitational darkening and also gives computational 
details. Section [3] describes the stellar and disk models chosen for this investigation. Re- 
sults are presented in Section H] where the consequences of rotation on the energy reaching 
the circumstellar disk (Section 14. ip and the corresponding changes in the disk temperature 
structure (Section 14. 2p are investigated. In Section [5l we examine the two principle ingre- 
dients of gravitational darkening (temperature variation with latitude and distortion of the 
stellar surface) separately in order to gauge their relative importance. We also evaluate the 
effectiveness of the simpler spherical gravitational darkening approximation in which the 
shape distortion produced by rapid rotation is ignored. Section |6] gives the conclusions, and 
in Appendix A, we discuss the subtle aspect of how varying the equatorial radius can change 
the inherent properties of our disk models. 



2. Theory of Gravitational Darkening 
2.1. The Central Star 



Gravitational darkening, first described by von Zeipel (1924), is due to a rotational, 
centrifugal term included in the classical gravitational potential. For a star rotating at a 
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fixed angular speed u, the potential, in spherical coordinates (r, 6, 0), is given by 

$fr,^) = - -wVsin"^, (11 



GM 1 2 2 . 2Z) 

UJ r sm t7, 

r 2 



where G is the gravitational constant, M is the stellar mass, and 6 is the stellar co-latitude 
[6 = 0° on the polar axis). The local gravitational acceleration is given by 



9{e) 



uo r sm 



in^ 



GM\ 



J r + (w^r sin6' cos^) 6 



(2) 



The angular speed at which the local value of gravity at the equator becomes zero defines 
the critical angular speed, Wcrit, and the corresponding critical rotation velocity, fcrit- The 
angular speeds and the equatorial velocities of stars are often expressed as fractions of these 
critical values, cufrac and Wfrac- The relevant equations are 



crit 



/SGM 



3r„ 



UJ 

feq 
^crit 



(3) 



where is the polar radius of the star. In our c alculations we have assumed that the polar 
radius remains constant following ICoUinsI ( 119661 ). 



By requiring the poter itial across th e surface to be constant, the radius as a function of 
is found to be (see again ICoUinsI Il966l ) 



r{e) 



-3rr 



Wfrac sin e 



( arccos (wfrac sin ff) + 47r 
cos^ 1 



(4) 



When the equator of the star becomes unbound at w = Wcrit, 



-- |) takes on its largest 

value, |rp. The Roche distortion of a stellar surface has direct interferometric confirmation. 



For example, IZhao et al.l (120101 ) present reconstructed images of a Cep and a Oph which 
clearly reveal rotational distortion consistent with the von Zeipel theory. The shape of the 
distorted surface affects the stellar radiation intercepted by the disk as the mid-latitudes 
are tipped towards the pole and away from the equatorial regions. In addition, rapid rota- 
tion increases the surface area and projected surface area of the sta r. The physical effec ts 
associated with near critical rotation have been recently reviewed by lMeynet et al.l (120101 ). 



The local effective temperature of the stellar atmosphere. Teg, and the local surface 
gravity, |^|, are related by von Zeipel's theorem. 



1/4 



(5) 
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where C^j is von Zeipel's constant. IVinicius et al.l (120071 ) find that the von Zeipel variation 
of Teg with stellar latitude is required to simultaneously fit several He I photospheric ab- 
sorption profiles from five rapidly rotating, early-type Be stars. Nevertheless, the power of 
1/4 appearing i n the standard gravity-darkening law has been observationally questioned by 
Monnier et al.l ( 120071 ) who fi nd a bette r fit to their interferometric observations of Altair with 
an exponent of 0.19±0.012. Ivan Belld ( l2010l ) reviews all current interferometric observations 
of rapidly rotating stars and finds that of six observed stars, four are well fit by the standard 
1/4 exponent and two (including Altair) are not. Despite these results, we will retain the 
standard exponent in all calculations to follow. 

The constant, C^j-, is defined using the assu mption that the tota l luminosity of the star 
remains fixed (for discussion of this point, see iLovekin et al.l 120061 ) . Therefore, while the 
equatorial temperature of a star decreases with rotation rate, the polar temperature must 
increase to maintain a constant luminosity. This effect is illustrated in Figure |T] which shows 
the change in temperature with increasing rotation for five stellar co-latitudes: 0°, SO'', GO'', 
70", and 90". In the absence of rotation, the stellar surface has a uniform temperature, 
but as the rotation speed increases the surface temperatures begin to diverge. Note that 
the temperatures of the middle co-latitudes, such as 60°, do not monotonically increase or 
decrease in temperature. For example as rotation increases from zero, the temperature at 
60° decreases until ~ O.SOwcrit where it reaches a minimum before increasing slightly with 
increased rotation. 

The constant can be found for a specific set of stellar parameters (r^, M and L), given 
a value of w, by integrating the local gravity over the stellar surface and then multiplying by 
the luminosity divided by the Stephan-Boltzmann constant. Once C^^ has been determined 
for one set of stellar parameters, the results can be rescaled for any other star following 
CoUinsI (119661) using 

' (6) 



M'L 



Here M , L , Sz are the mass, luminosity and von Zeipel's constant of a particular star, 
and M, L, & C^^ are for another star with an identical rotation rate. 

The local gravity and surface orientation can be established from Equation [2l and the 
temperature from Equation [51 This allows the properties of the local stellar atmosphere, g 
and Teff, to be defined for each point on the surface of the star. 
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Fig. 1. — Change in the local effective temperature with increasing ffrac (lower axis) for a 
B2V stellar model (see Tabled]). The upper axis shows the corresponding Ufrac- The darkest 
solid line gives the polar temperature [6 = 0°), while the lightest grey dash-dot line shows 
the equatorial temperature {6 = 90°). Various intermediate stellar latitudes are also shown. 
The right vertical axis shows the standard main sequence s pectral type (with no rotation) 
for the same surface temperature as the left axis (ICoxll2000l ). 
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2.2. The Circumstellar Disk 

We now turn to the Be star circumstellar disk. From a particular vantage point within 
the disk, labelled s, defined using the cylindrical coordinates, {sr, Sz, only a portion of 
the star is visible. For a section of the stellar surface specified by r{6, (p) to be visible from 
the vantage point, the dot product between the unit vector pointing from the stellar surface 
to the vantage point, h^p = tt^, and the local surface normal, -n-surf = -nr, must be greater 
then zero. The boundary of this region where ngurf ■ n^p = is important because it allows 
the lines of sight between the star and the vantage point to be efficiently chosen. The upper 
and lower edges of the visible region correspond to the maximum and minimum values of 6 
such that risurf ■ "^vp > 0. These are found by solving 

8 

T^^I^^c-^sr sin 9 - rg - SRsinO + cos 6^ = (7) 
27 

for it's two solutions, di and 92- For values of 6 between di and the visible portion of the 
surface lies between the bounds on given by 

, 'r(6') sin^ 6'(1 - Kg) + r(6') cos^ 6* - cos 61 

'edge = s$ ± arccos 



(s/j sin 9 — KgSji sin 6) 
with ^ 

Given this geometry, we can now take into account the temperature changes and the 
radial distortion of the star to determine the photoionizing radiation field reaching each 
point in the disk. The effects of gravitational darkening are illustrated in Figure [2l This 
figure shows a star rotating at 80% of Vcrit with its shape distortion and associated surface 
temperature variations. Figure |2] also shows the visible region, 9i, 62 and 0edge, for a vantage 
point at a radius of = 2 and a height above the disk of S2 = 1.2 r^. The disk is assumed 
to be axisymmetric so any $ dependence is not included. This figure demonstrates how the 
visible region changes from a simple spherical cap for a non-rotating star to an elongated 
region. 

The circumstellar gas is exposed to photoionizing radiation from the various latitudes 
of the star and this radiation changes with increasing rotation. Different locations within 
the disk have different regions of the photosphere within their field of view. The polar 
temperature increases with rotation in contrast to the decrease in the equatorial temperature 
(see Figured]). The drop in the equatorial temperature occurs at a faster rate than the rise 
in the polar temperature. This means the effects of the hot stellar pole become significant 
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at speeds above 80% of fcrit- In the equatorial plane up to a distance of 4rp, stellar latitudes 
above 25° are not visible for a circular star. For a critically rotating star this increases 
to ~ 30°. As the gas density in the disk is assumed to fall off with height above the 
equatorial plane, the light reaching upper parts of the disk from the cooler stellar equatorial 
region passes through a greater optical depth than that from the hotter pole. Thus the 
effects of gravitational darkening can vary in strength depending on location within the disk 
and can have very different effects on the temperature depending on what part of the disk 
is considered. The effect of gravitational darkening cannot be considered to be simply a 
reduction of the effective temperature of the star. 



3. Calculations 



A gravitationally darkened version of Bedisk (jSigut fc Joned 120071 ). incorporating the 
theory of Section [21 was created and run for the stellar parameters given in Table [I] and 



the rotation rates given in Tabled These stellar parameters (adopted from Coxll2000l ) were 
chosen to include a model from four of the five bins used by ICranmerl (120051 ) to analyze the 
effects of spectral type on Be star rotational statistics. For reasons discussed below, we have 
not considere d models f rom t he coolest bin (the bin with characteristic spectral type B8) 
considered by ICranmerl (120051 ). Bedisk solves the statistical equilibrium equations for the 
atomic level population equations and then enforces radiative equilibrium at each point of 
the computational grid ( ISigut fc Joned 120071 ). The density structure of the disk is assumed 
to be of the form 



p{R, z) = po \ — \ e 



~R 



(10) 



Here R^. is the stellar radius (see next paragraph). In the second equation for the disk scale 
height, H{R), p is the mean-molecular weight of the gas, and Tiso is an assumed isothermal 
temperature used for the sole purpose of setting the density scale height. The models are 
constructed with the density parameters set to n = 3 and po = 5 x 10~^^ g/cm^. These 
density parameters are kept constant. The form of the disk density given by Equation [10] 
follows from the assumption of a radial power-law drop in the equatorial plane [z = 0) 
coupled with the assumption that the disk is in vertical hydrostatic equilibrium set by the 
^-component of the star's gravitational acceleration. Note that Equation [TU] results in a 
flaring disk with H oc i?^/^. 
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Fig. 2. — The temperature structure of a B2V star rotating at ffrac = 0.8 viewed with an 
inchnation of 45° to the polar axis (marked by a small dot in uppermost oval). The bold, 
solid line is the equator. The local temperatures shown are determined by Equation (j5]). 
The dashed line is the limit of the region visible from a vantage point located within the 
disk at a radius of sr = 2rp and a height of = 1-2 r^, as defined by Equation ([8]). The 
square and diamond symbols mark the location of 6i and 62 respectively which are found by 
solving Equation (JTj). The direct line of sight from sr = 2rp and Sz = 1.2 rp is indicated by 
the symbol 0. 
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If the equatorial radius expands with rotation, this can change the disk's density struc- 
ture if the appearing in Equation [10] is associated with the star's equatorial radius. This is 
undesirable because we wish to investigate changes in the disk temperatures resulting solely 
from the changing photoionizing radiation field seen by the disk and not from an unintended 
change to the underlying density structure of the disk. Unfortunately, there is no ideal solu- 
tion to this complication and several approaches are outlined in Appendix A. In the notation 
of Appendix A, all of the calculations were computed with an unchanging grid in which the 
inner edge of the disk is set to 3/2 rp. This approach preserves the disk density structure 
and is best for temperature comparisons. A physical grid, in which the density structure 
of the non-rotating star is shifted outward as i?^, increases, is more useful when calculating 
observables, and this grid will be used in a subsequent paper on predicted Ha line profiles 
and infrared excesses. 



Table 1: Adopted stellar parameters 



Spectral 


Mass 


Polar Radius 


Luminosity 






3 

^crit '^crit 


Type 


(Mo) 






(s-1) 




(kms^"*^) 


BOV 


17.5 


7.40 


3.98x10^ 


7.10x10" 


-5 


548 


B2V 


9.11 


5.33 


4.76 xlO^ 


8.38x10" 


-5 


466 


B3V 


7.60 


4.80 


2.58 xlO^ 


8.95x10" 


-5 


449 


B5V 


5.90 


3.90 


7.28 xlO^ 


1.08x10" 


-4 


439 



Notes: Stellar parameters adopted from ICoxl (120001 ) 



Table 2: Adopted rotation rates and example equatorial speeds. 



^frac 


<^frac 


req/rp 


^^eq BO 


Vcq B3 


Vcq B5 








(kms~^) 


(kms^^) 


(kms^^) 


0.000 


0.000 


1 











0.001 


0.002 


1.00 


0.55 


0.45 


0.40 


0.200 


0.296 


1.01 


110 


90 


80 


0.400 


0.568 


1.06 


219 


180 


160 


0.600 


0.792 


1.14 


329 


269 


241 


0.800 


0.944 


1.27 


438 


359 


321 


0.950 


0.996 


1.43 


521 


427 


381 


0.990 


0.999 


1.49 


544 


445 


397 
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The photoionizing radiation field incident on the disk was calculated as follows: first, 
von Zeipel's c onstant is found by using Equation [6]to rescale an interpolated value of from 



CoUinsI (Il966[ ) to match the stellar parameters. Next, for each point in the disk, the visible 
region of the stellar surface is determined, and a grid is constructed across this surface. At 
each point in this surface grid, the local temperature, local gravity, and the viewi ng angle 



are fo und. Finally, the specific intensity, ly (Teff,g), for each point is scaled from the iKurucz 



( I1993I ) model atmosphere closest to the effective temperature and gravity at each location, 



(Tbest, 5'best), to the required temperature using a blackbody function, 

If {TeB, g) = (^best, 5'bcst) p "/^''^l • (H) 

Ubcstj 

Here Tcfj is the local effective temperature of the stellar surface, and Tbest and (?best are 
the parameters of the closest matching Kurucz model. The available Kurucz grid of model 
atmosphere intensities was interpolated to a finer grid of spacing of ATefj = 200K and 
Alog((y') = 0.10. Equation [TT] was used to correct the intensity only in these very narrow 
temperature bins. 

Most of the calculations performed in this work implement the above description of 
gravitational darkening and are referred to as full gravitational darkening (FGD) calculations. 
All of the models presented in Section H] are FGD calculations. However in Section [5l two 
simplified approaches are described which isolate one of the two main effects of gravitational 
darkening: spherical gravitational darkening (SOD), accounting only for the temperature 
change across the stellar surface, and pure shape distortion (PSD), accounting only for the 
changing the shape of the star. 



4. Results for Gravitationally-Darkened Be Star Disks 

4.1. Energy reaching the Circumstellar Disk 

The stellar radiation field is modified by gravitational darkening, and this directly af- 
fects the disk temperatures through the photoionization rates. Figure [3] shows the total 
photoionizing mean intensity at each location in the computation grid for a B2 star both 
with and without a disk, and with and without rotation. The mean intensity ratios between 
the various cases are also shown. Consider the left-hand panels A and B in which there is no 
disk material. In the rotating case (ffrac = 0.99), less radiation from the poles reaches the 
equatorial regions of the disk, and this results in the hour glass shape in the mean intensity 
of panel B. In the bottom-left panel which gives the ratio of the rotating to non-rotating 
case, the stellar radiation field reaching the disk has been reduced by approximately 50% 
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compared to the non-rotating star. However, for regions close to the star and near the equa- 
torial plane, there can be reductions as great as 90%. Now consider the right-hand panels 
C and D in which the opacity and emissivity of the disk gas have been included. Again 
the cases of rotation at i^frac = 0.99 and no rotation are compared. Most of the radiation 
reaching various positions within the disk is re-processed by the circumstellar gas between a 
particular location and the star. In general, the total irradiance reaching the shielded parts 
of the disk is reduced down with gravitational darkening by about 50%. However, there are 
locations within the mid-regions of the inner disk where the energy actually increases with 
rotation (see Figure [3]) due to changes in the disk's opacity. 



4.2. Temperatures in the Circumstellar Disk 

Any change in the energy received from the star produces a change in the radiative 
equilibrium temperature. The changes produced by increasing rotation on the general tem- 
perature structure of the disk are illustrated in Figures H] through [7] for the spectral types 
given in Table [H Shown in all figures are four disk temperature diagnostics: the maximum 
temperature, the minimum temperature, the density-weighted temperature, defined as 

^P = jr- [ TiR,z)piR,z)dV , (12) 

and the volume-averaged temperature, defined as 

Tv = -^ / T{R,z)dV. (13) 

Vdisk J 

In order to avoid numerical effects, the maximum and minimum temperatures are averages 
over the twenty hottest and twenty coolest disk locations, respectively. In Figures H] through 
[7l all of these temperature diagnostics are plotted as a function of ffrac and Wfrac- Because 
their behaviour depends somewhat on spectral type, the results for each model are discussed 
separately. 

For the BO model, seen in Figure HI both the density and volumed-averaged disk tem- 
peratures decrease steadily with increasing rotation. Small declines in the density-weighted, 
maximum, and minimum disk temperatures are apparent even at ffrac ~ 0.4. By critical 
rotation, the density-weighted temperature has fallen to just under 11, 000 K, a decline of 
2500 K compared to the non-rotating case. The decline is not as large in the volumed-average 
temperature, and this is likely due to the influence of the hotter stellar pole that develops for 
increased rotation (see Figured]); this is also reflected in the maximum and minimum disk 
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Vacuum Disk 




5 10 15 20 5 10 15 20 

R/r R/r 
p p 




Fig. 3. — Left panels A and B show the log of the mean intensity at each location in the 
computational grid without any disk material present for a star with ffrac = (panel A) 
and ffrac = 0.99 (panel B). The ratio of the mean intensity for the rotating case with the 
non-rotating case is given in the bottom left panel (labelled B/A). The right panels C and 
D show the mean intensity at each location in the disk including the opacity effects of the 
disk material for a star with ffrac = (panel D) and ffrac = 0.99 (panel E). The ratio of the 
shielded rotating case with the non-rotating case is shown in bottom right panel (labelled 
D/C). 
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Fig. 4. — Change in disk temperature diagnostics with increasing rotation for the BO modeL 
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Fig. 5. — Same as Figure H] for a B2 star. 
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Fig. 6. — Same as Figure H] for a B3 star. 
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Fig. 7. — Same as Figure H] for a B5 star. 
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temperatures: while the minimum temperature has a steep dechne, the maximum temper- 
ature is much flatter and actually reaches a shallow minimum at an intermediate rotation 
rate. 

For the B2 model, shown in Figure [5], all temperatures initially decline slowly. The 
density- weighted disk temperature falls from over 9000 K in the non-rotating case to under 
8000 K by critical rotation. However, unlike the previous case, the volume-averaged temper- 
ature reaches a minimum near ffrac ~ 0.8 and then begins to increase for extremely rapid 
rotation. The lack of an increase in the density-weighted temperature indicates that the 
additional heating occurs mainly in the upper disk, where there is less gas, and the heating 
is likely due to the effect of the hot stellar pole. We also note that minimum temperature in 
the disk, ^ 6000 K, is approximately constant with rotation, with only a small decline for 

^^frac > 0.7. 

Spectral type B3 is shown in Figure [61 Its behaviour is very similar to the B2 case 
above, although the temperature rise in the volume-weighted average temperature above 
"Wfrac > 0.8 is not as large as in the B2 case. Again the minimum disk temperature is around 
6000 K and is not strongly affected by rotation. For B5, shown in Figure [TJ both the density 
and volume-averaged temperatures show a similar decline with rotation and there is no rise 
in the volume-averaged temperature above ffrac > 0.9. By this spectral type, the minimum 
disk temperature has fallen to ~ 5000 K and shows no clear trend with rotation. The low 
temperatures reached in this model have implications for later spectral types that we shall 
now discuss. 

As noted in the introduction to Section [3l we have not co nsidered a model of spectral 



type B8 from the VL (very late) spectral bin of ICranmen (120051 ). Our test calculations have 
indicated that very cool disk temperatures are reached in this model, particularly for near 
critical rotation, and these may fall outside of the domain of applicability of the current 
version of BEDISK. In particularly while the global density-weighted disk temperature still 
exceeds 6000 K for a B8 model, there is a significant volume of the disk which falls well 
below 5000 K by t>frac = 0.99. While our calculations include abundant metals with low 
ionization potentials (such as Mg, Ca, and Fe) that provide sources of free electrons at 
low temperatures, molecule formation is not included. In addition for such cool disks, the 
treatment of the diffuse radiation field generated within the disk may require a more careful 
treatment. For these reaso ns, we have cho sen not to include B8 models in this work. We 



note that in the analysis of ICranmed (120051 ). the four bins we do account for include the full 
range of inferred threshold rotation rates, from ffj.ac~0.6 for early spectral types B0-B2 to 
t)fi,ac~1.0 for the late spectral type B5. 

The effect that near critical rotation (ffrac = 0.99) has on the global, density-weighted 
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disk temperature is summarized in Table [31 Compared to non-rotating models, near critical 
rotating models are between 15 and 20% cooler with the largest difference occurring for the 
earliest spectral type considered (BO). 

Next we discuss how the detailed temperature structure of a Be star disk, T{R,z), 
changes with rotation. In general. Be star circumstellar disks are highly non-uniform in 
temperature as disks of sufficient density develop a cool, inn er zone surrounded by a larger 



warm region (IMillar fc Marlborough! Il998l : ISigut et al.ll2009l ). Between these two regions is 



the hottest part of the disk, forming a narrow sheath. Figures [Hand [HI show the effect of rota- 
tion on the temperature of disks surrounding stars of spectral types BO and B5 respectively. 
In each plot, temperature distributions for two models are shown, one with t>frac = 0.0 and 
the other with ffrac = 0.95. Also shown is the temperature difference between the rotating 
and non-rotating models and a histogram of these differences. The value of ffrac = 0.95 was 
chosen for the rotating model because this is approximately the rotation rate required for a 
velocity perturbation on the order of the local sound speed in the star's photosphere to be 
sufficient to feed material into a rotationally supported disk. 

For the BO model, shown in Figure [HI the most dramatic changes in temperature occur 
because of the expansion of the cool, inner zone to larger radii with increased rotation. In 
the non-rotating case, this cool region does not extend past ~ 7i?* while in the ffrac = 0.95 
model, it extends well past 10 i?*. This results in some computational grid points that 
are more than 4, 000 K cooler in the rotating model (compared to the 2500 K difference in 
the density-weighted average) as illustrated by the histogram of temperature differences in 
Figure [H] However, there are also temperature changes in the hot sheaths above and below 
the inner cool zone close to the star. For the ffrac = 0.95 model, these hot sheaths occur only 
beyond ~ 3 i?* whereas in the non-rotating case, these sheaths extend to the inner edge of 
the disk. 

Table 3. Change in density-weighted average disk temperature for near-critical rotation. 



Spectral Type Tp{K) % Change 

^^frac = 0.00 t;frac = 0.99 



BO 13280 10780 -19 

B2 9190 7870 -14 

B3 8690 7460 -14 

B5 7660 6760 -12 
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Fig. 8. — Temperatures for both a non- rotating BO star and its circumstellar disk (upper 
left) and a BO star rotating at ffrac = 0.95 with an identical disk (lower left). In each of these 
panels, the colour bar on the left is for the star and the colour bar on the right, for the disk. 
The lower right panel shows the temperature differences in the disk. Positive differences 
mean that the non-rotating star is hotter. The upper-right panel shows a histogram of the 
temperature differences. 
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Fig. 9.— 



Same as Figure [8] for a B5 star. 
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Despite these changes, the temperature in the equatorial plane of the cool, inner zone 
does not drop dramatically with rotation, likely due to the high optical depths along the rays 
back to the star. This explains why the density-weighted average temperature of the disk 
shows only a modest decrease compared to the non-rotating case (see Figure H]). Interestingly, 
the rotating model is actually hotter than the non-rotating model in a narrow region a few 
scale heights above the equatorial plane for R > lOi?*, although the number of computational 
grid points involved is quite small. The optically thin gas above and below the plane of 
the disk shows only a very small decrease in temperature and thus the volume-averaged 
temperature of Figure H] is essentially unaffected by rapid rotation. 

Figure [8] clearly illustrates that the temperature distribution surrounding a rotating 
star may differ dramatically from the non-rotating case, even when global measures of disk 
temperature such as Tp and Ty are similar. From the point of view of computing the strength 
of emission lines or the infrared excess, the increased extent of the inner cool zone and the 
changes to the temperature of the hot sheaths above and below this zone have the potential 
to produce large changes in these observational diagnostics. Such changes will be discussed 
in a subsequent paper. 

Figure [9] shows a similar comparison between a non-rotating model and one rotating 
at ffrac = 0.95 for spectral type B5. As in the BO case, the main effect is the extension of 
the inner, cool zone to larger radii. However the B5 model also has a significant extension 
of this cool zone to larger distances above (and below) the plane of the disk. As these 
regions occupy a significant volume, the volume-weighted disk temperature is significantly 
decreased (unlike the previous case of a BO star). As shown in Figure [TJ both the density 
and volumed-weighted disk temperatures show a similar drop with rotation. 

While most disk locations are, as before, cooler, there is a much larger fraction of disk 
locations that are hotter in the B5 rotating model than in the previous BO case. Most of 
these hotter regions are associated with a cool zone at large radii, R > 10 R^, far above 
(and below) the equatorial plane of the disk. In these regions, the lines of sight back to the 
central star have significant optical depths whereas the optical depths (in the 2;-direction) 
to the nearest disk edge have become less than unity. Hence heating via photoionization 
is suppressed whereas cooling via escaping line radiation becomes effective and there is a 
reduction in temperature due to this cooling. This low temperature zone disappears for 
higher z because the optical depths back to the central star eventually drop and the heating 
rate increases. The exact location and detailed form of the zones is influenced by rotation 
and some locations with significantly hotter temperatures are predicted. It should be kept 
in mind that at large radii and far above or below the plane of the disk, the gas densities are 
very low, and the detailed structure of these cool zones will have little effect on observable 
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diagnostics such as emission line strengths, polarization, or infrared excess. 

While the detailed T{R, z) comparisons of the previous paragraphs are instructive, it 
is difficult to see general trends over a wide range of rotation rotates. For this reason, it 
is useful to have temperature diagnostics that are intermediate between the disk-averaged, 
global measures of Figures H] to [7] and the detailed temperature distributions of Figures [8] and 
m To this end, we have performed both vertical {z) and radial (i?) temperature averages 
over the disks of three stellar models (BO, B3, and B5) and considered how these averages 
are affected by a wide range of rotation rates. 

The vertical temperature averages (yielding an average disk temperature as a function 
of distance from the star) were computed via 

, rf™'''' TiR, z)p(R, z)dz , ^ 

Jo P{R^ dz 

and the change in vertically-averaged, density-weighted temperature as a function of radius 
is shown in Figures fTOl through [T2] for three stellar models, BO, B3, and B5. Five rotation 
rates, ffrac = 0.00, 0.40, 0.80, 0.90, 0.95 and 0.99, were considered. The behaviour of the B2 
model of Table [H not shown for brevity, is similar to the B3 case. In all cases considered, 
the ffrac = 0.40 models show vertically averaged temperatures that are very close to the 
non-rotating case while significant differences develop by ffrac = 0.80. 

These vertical averages show that the region next to the central star is hottest, followed 
by a rapid decrease towards a minimum temperature near ^ 4 R^, (for BO) to ~ 6 i?* (for B5). 
Further out, the temperature increases again towards the outer region of the disk. In the BO 
model, the temperature then plateaus, whereas in the B3 and B5 models the temperature 
is still increasing by i? ~ 20 i?*. From these figures, it can be seen that increasing rotation 
tends to: (1) decrease the average temperature at all radii, (2) move the location of the 
temperature minimum in the cool zone inward (although not in the BO model), and (3) 
broaden the extent of the cool, inner zone. Rotation also tends to fiatten the temperature 
increase occurring in the outer regions of the disk. 

Turning to the vertical {z) structure of the disks, we define a density-weighted, radially- 
averaged temperature (yielding an average temperature as a function of height above the 
equatorial plane) as 

J^''--TiR,u)p{R,u)A{R,u)dR 

Tiu) = — p , (15) 

J^'-^'' p{R,u)A{R,u)dR 

where u = z/H{R) and H{R) is the vertical scale height given by Equation [101 The vertical 
variable, z, was rescaled by H{R) because lines of constant u better define a radial average 
due to the increase in the disk scale height with R (see Equation [TOj) . The function A{R,u) 
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Fig. 13. — Variation of the radially averaged temperature profile in the vertical direction for 
the BO model. Here the gas scale height is the variable u defined hj u = z/H{R). 
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arises because Equation [15] is effectively a line integral through the computational grid. At 
each Ri in the discrete sum used to compute this integral, z must be found from Ri and the 
required u. As different z occur in the sum and the spacing of the z-grid increases with R, 
we have weighted each z point by the width of the strip it effectively represents as based on 
the local grid spacing; this is the origin of the A{R, u) function. Such considerations do not 
occur for Equation [H] because this integral is performed over all Zj for a fixed Ri. 

The changes in the vertical temperature structure with rotation are shown in Figures [T5] 
through [TBI for spectral types BO, B3, and B5, respectively. Again, the B2 model of Tabled] 
is similar to B3 but is not shown for brevity. Note that for each spectral type, the set of 
weights A{R, u) used in Equation [15] is always the same as the computational grid is fixed 
and independent of ffrac (i-e. we have chosen to use the unchanging grid of Appendix A). 

In general with increasing vertical height (i.e. perpendicular to the equatorial plan), the 
temperature is coolest in the plane of the disk [z = 0) and then rises to a maximum near 
u~ 2. This maximum corresponds to the hot sheath seen in the plots of the temperature 
distributions in these disks (for example, see Figure [H]). At larger heights, the temperature 
drops below the mid- plane maximum and is nearly constant. With increasing rotation, 
several trends are seen: (1) rotation reduces the temperature at all scale heights, (2) rotation 
increases the vertical extent of the cool, inner zone, (3) rotation decreases the temperature 
maximum in the mid-plane, and (4) rotation shifts the maximum temperature to larger scale 
heights, although this shift is not large. 

Nevertheless, there are exceptions to these general trends, particularly in the tempera- 
tures at large u. In the BO model (Figure [T3|) . the temperature at the upper edge of the disk 
decreases monotonically with rotation. However, for the B3 model, the temperature at the 
upper edge of the disk is actually largest for the greatest rotation rate. The location of the 
hot sheath (the temperature maximum in u) moves outward in u for later spectral types; it 
occurs well below u = 2 for the BO model at all rotation rates, whereas it is above u = 2 for 
the two later spectral types, and it moves higher with increased rotation. Hence in the later 
spectral types, the temperature at the upper edge of the disk is affected by the presence of 
the hot sheath. 

Finally, turning back to global measures of disk temperature. Figure [T6] summarizes the 
change in the density-weighted average disk temperature with rotation for all four spectral 
types. In general, the non-rotating models define the upper envelope of the temperatures 
while the extreme rotators (ffrac > 0.95) define the lower temperature envelope. In con- 
structing this figure, any change in the overall spectral type of the star due to increasing 
rotation has been ignored. 
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5. Approximate Treatments of Gravitational Darkening 

In this section, we assess the effectiveness of partial or approximate treatments of grav- 
itational darkening. In one approach, spherical gravitational darkening (SGD), we have 
treated the star as a sphere and implemented only the variation of the photospheric tem- 
perature across the stellar surface. As this procedure is very simple to implement, it is of 
interest to see how the predicted disk temperature structure compares to the full treatment. 
In another approach, pure shape distortion (PSD), we include only the rotational distortion 
of the central star, and not the temperature variation, to understand the effects of geome- 
try alone. Global density and volume-weighted disk temperatures for both approximations 
(SGD and PSD), and the full treatment (FGD), are shown in Figure [T71 for the B2 stellar 
model. 

The differences between FGD and SGD are subtle for Vf^ac < 0.7, with only small 
differences in the density-weighted and volume-weighted average temperatures. For the 
largest rotational rates, ffrac > 0.95, SGD predicts density-weighted temperatures that level 
off for higher rotation rates while the full treatment gives temperatures that continue to 
drop. The difference amounts to ~ -1-500 K for the SGD model at Wfrac = 0.99 in the density- 
weighted average temperature. In the volume-averaged temperature, the FGD temperature 
actually falls below the SGD prediction in the vicinity of ffrac ~ 0.80, but then increases 
above the SGD prediction for larger rotation rates. By ffrac = 0.99, the difference is again 
about 500 K. 

To further illuminate this result, these calculations were repeated with PSD. As evident 
from Figure [T71 PSD results in essentially no change to the average disk temperatures for 
■^frac < 0.7 and only very small changes for faster rotation rates. This, of course, does not 
mean that the effect of the rotational distortion of the stellar surface is unimportant; indeed 
the difference between the FGD treatment and SGD noted above is the neglect of the stellar 
distortion. The PSD treatment simply illustrates that the distortion of the stellar alone 
produces almost no change in global disk temperatures. 

However global temperature diagnostics tell only part of the story. To illustrate how 
the temperature structure of the disk is reproduced by the SGD approximation. Figure [T8] 
compares the vertically-averaged (Eq. [H]) and radially-averaged (Eq. [15]) disk temperatures 
for a B2 stellar model computed with three rotation rates, ffrac = 0, 0.80, and 0.95. As can 
be seen from the Figure, there are significant differences between all the rotating models 
(either FGD or SGD) and the non-rotating model, and all of the previously discussed effects 
can be seen. Interestingly both SGD profiles for ffrac = 0.8 are close to the FGD prediction. 
However by ffrac = 0.95, the SGD temperature profiles remain close to the ffrac = 0.8 SGD 
predictions and do not follow the trend of the full (FGD) treatment; this is particularly 
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noticeable in the radially-averaged profile (bottom panel) where the (average) location of 
the hot sheath is very poorly predicted by the SGD model. We conclude that SGD is an 
acceptable approximation to the temperature structure only for i^f^ac < 0.8, and that for 
higher rotation rates, the full treatment is required. 

6. Conclusions 

Gravitational darkening produces noticeable changes in the temperature structure of Be 
star circumstellar disks. Rotation causes the stellar surface temperature to diverge from a 
single value into a range of temperatures, and this changes the photoionizing radiation field 
incident on the disk. Three measures of disk temperature were used to quantify the effect of 
rapid rotation: detailed temperature distributions, T{R,z), vertically and radially-averaged 
disk temperatures, and global density or volumed-weighted average disk temperatures that 
reduced the complex temperature variations of the disk to a single, average temperature. 
Among these choices, the vertically and radially-averaged temperatures are perhaps the 
most useful for illustrating the effect of rapid rotation over a wide range of rotation rates 
and spectral types. 

Rapid rotation of the central star and the accompanying gravitational darkening cause 
most of the disk gas to systematically decrease in temperature, and this is reflected in lower 
density-weighted, average disk temperatures for increasing rotation rates. However volume- 
averaged disk temperatures indicate that upper regions of the disks can experience additional 
heating at rotational speeds above 80% of the critical rate. Such changes cannot be produced 
by simply lowering the effective temperature of the star with rotation. 

All spectral types considered, B0-B5, show a decline in the global, density- weighted 
averaged disk temperature of ~ 15-20% for near critical rotation. The decline is largest for 
the earliest spectral type considered, BO. Changes in the global, volume-weighted average 
disk temperature and the minimum and maximum temperatures in the disk also occur, 
although results are dependent on spectral type. The volume-weighted average temperature 
shows the effect of the hot stellar pole (that develops for rapid rotation) at intermediate 
spectral types B2 - B3 (see Figures |5] and [6]) . 

The most important change to the temperature structure of the disks is the expansion 
of the inner cool, zone close to the central star. This zone develops in all disks of sufficient 
density, and this cool region generally expands in both radius and height as the rotation rate 
is increased. Also important are the hot sheaths above and below the cool, inner zone. These 
generally decrease in temperature and move to larger scale heights with increased rotation. 
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Fig. 16. — Density-weighted temperature average plotted against spectral type for six dif- 
ferent rotation rates. 




Fig. 17. — The global volumed-averaged and density-averaged disk temperature as a function 
of ffiac for the B2 model. Three different treatments of gravitational darkening considered: 
FGD, SGD, and PSD. 
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Fig. 18. — Vertically-averaged (top panel) and radially-averaged (bottom panel) temper- 
ature profiles (see Equations [H] and [T5|l for a B2 stellar model that compare the SGD 
approximation with the full treatment (FGD). Models rotating at Vfrnc = (solid lines), 0.80 
(dot-dashed lines) and 0.95 (dashed lines) are shown. 
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It is these complex changes that will most affect physical observables such as the Ha emission 
line strength or the infrared excess. These issues will be explored in a subsequent paper. 

The temperature effects of gravitational darkening on global measures of Be star cir- 
cumstellar disks can be adequately approximated by spherical gravitational darkening (SGD, 
which ignores the distortion of the stellar surface) for rotation rates less than 80% of critical. 
For higher rotation rates, the distortion of the stellar surface must be included to accurately 
compute the disk's global average temperature. Comparisons of the vertically-averaged and 
radially-averaged disk temperature profiles for models with different rotation rates suggests 
that the general temperature structure of the circumstcllar disk can be reasonably predicted 
using the SGD approximation for rotation rates of also less than 80% of critical. 

While the effect of rapid rotation of the central star can produce complex effects in the 
temperature structure of Be star circumstcllar disks, temperature is, of course, not directly 
accessible to observers. In a subsequent paper, we will quantify the effect of rapid stellar 
rotation on the classic diagnostics of circumstcllar material, namely emission line strengths 
and profiles and the predicted near-IR excess in the system's spectral energy distribution. 
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A. Construction of the Computational Grid 

The primary goal of this paper is to explore changes in the thermal structure of a 
circumstellar disk caused by changes in the stellar flux due to rotation. Ideally the only 
difference between different models should be the rotation speed of the star. Rapid rotation 
has two effects on the central star; the dependence of the photospheric temperature on stellar 
latitude, and the Roche distortion of the stellar surface. If the latter effect (distortion) did 
not occur, creating a series of comparison models would be a simple matter of swapping out 
the central star and examining the corresponding changes in the thermal structure of the 
disk. Indeed, this is what occurs in the SGD approximation. Nevertheless, the stellar surface 
is distorted by rotation, and this geometrical change can be propagated into the physical 
properties of the disk. The simple fact that the distance between the stellar pole and equator 
changes with rotation is the central issue. For a given point within the disk it is possible to 
preserve the distance to one, but not both of these locations. Given this, there are several 
different ways one could attempt to mitigate the effects, and these approaches are discussed 
in this Appendix. All of the various choices are illustrated in Figure ^\ which shows in the 
left hand panels, a non-rotating B2 star and its associated {R, z) grid and in the right hand 
panels, a B2 star rotating at ffrac = 0.99 and its associated (i?, z) grid. In each panel, the 
(i?, z) grid points in the circumstellar disk are shown as small dots. 

In the most straight forward description, the grid locations (i?, z) within the disk are 
specifled in terms of the equatorial radius. In the absence of rotation, this is same as the 
polar radius. However when r^q increases due to rotational distortion, the computation grid 
is stretched in both R and z (see Eq. [10]). This increases the volume and mass of the disk. 
The computational grid points are also moved further from both the stellar pole, the stellar 
equator, and each other, as rotation increases. Hence with grids constructed in this way, any 
comparison between global temperature averages and total emission for increasing rotation 
must be considered carefully. Because of the simplicity of this approach, it will be referred 
to it as the naive grid. It is illustrated in panels A and B of Figure [191 

A simple alternative would be to deflne the grid spacing and density structure to be 
that of an non-rotating star, keep the grid flxed, and then shift the whole grid and density 
structure outward as req increases. This preserves the distances of all grid points from the 
stellar equator and between each other. While the cross-sectional area and density structure 
remains unchanged, both the inner and outer radii of the disk are systematically increased 
causing the volume and mass of the disk to also increase with rotation. This effect can 
be noticed when the disk emission measure is calculated. The distance to the pole of the 
star also increases systematically with rotation, and this could bias any global temperature 
averages to lower temperatures. This approach will be referred to as the physical grid and 
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is illustrated in panels C and D of Figure [191 

An alternative to maintaining the distance between the grid points and the equator is 
to maintain the distance between grid points and the polar axis. Once again the grid spacing 
and density structure is defined to be that of an non-rotating star, but instead of shifting 
it outward as req increases we simply remove sections of the grid that would be swallowed 
by the star and leave the rest of the disk unchanged. Unfortunately, the inner region which 
gets swallowed is also densest region of the disk. More of this region is removed as rotation 
increases, and this systematically reduces all optical depths within the disk which in turn 
increases the stellar fiux reaching the outer regions. This approach will be referred to as the 
swallowed grid and is illustrated in panels E and F of Figure [191 

One potential solution to the swallowed grid is to start the disk at R = (3/2) which 
allows the central star to swell to the inner disk boundary by critical rotation; the disk 
remains truly unchanged. This preserves disk volume, mass and density structure, and 
keeps the distances between each grid point and the stellar pole, the hottest part of the star 
unchanged. However, the distances to the stellar equator still systematically change. We call 
this the unchanging grid and it is illustrated in panels G and H of Figure [THl This method, 
unfortunately, makes important regions of the disk empty that were previously full and while 
comparisons between runs are unbiased any comparison to previous work becomes difficult. 
The dense inner region of the disk is the source of the IR lines and a significant amount of 
the mass in this region is no longer included. 

In conclusion, there are four ways to define a grid surrounding an expanding star: a 
naive approach where everything increases, a physical grid which shifts the disk, a swallowed 
grid in which mass is lost from the disk; and an unchanging grid in which the star closes an 
inner gap. The results included in this paper were computed with the unchanging grid to 
allow the most straightforward temperature comparisons between the difference models. 
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Fig. 19. — Each panel shows a cross section of the disk with the [R, z) grid points shown as 
small dots. The spherical star is not rotating, and the distorted star is rotating at Wfrac = 0.99. 
Panels A and B depict the naive grid, panels C and D, the physical grid, panels E and F, 
the swallowed grid, and G and H, the unchanging grid. The colours represent the log of the 
disk density. 



